% Matlab code for the numerical analysis of the paper "A common thread linking the design of 
% guarantee and non-escalating payments of public annuities" by Lau and Zhang (2023)
% Jan 14, 2023
clear all

% ------------------------ common paratmers ------------------------------% 
r = 0.20;
rho = 0.15;
phi = 2;  % CRRA parameter
zeta = 0.9;
m = 3.15;
psi = 0;  % correlation coefficient

wlb = 2.5; 
wub = 3.5;      
mu_w = 3;
eps2_w = 0.15^2;  
acr2 = 0.01;     % approximzation accurancy for wage rate distribution;
w = wlb:acr2:wub;

    % -------------------------  two cases ------------------------------
for i=1:2
    if i == 1    % limited heterogeneity case 
        thetalb = 0.2; % survival probability, lower bar
        thetaub = 0.3; % survival probability, upper bar  
        mu_theta = 0.25;
        eps2_theta = (0.0141)^2;
        acr1 = 0.0001;  % approximzation accurancy for health distribution;
        display(' ');                
        display(['Limited heterogeneity case: theta=[',num2str(thetalb),',',num2str(thetaub),'], psi=',num2str(psi),', phi=',num2str(phi),', zeta=',num2str(zeta),'.']);        
    elseif i == 2   % substential heterogeneity case 
        thetalb = 0.5; 
        thetaub = 0.99;  
        mu_theta = 0.7;
        eps2_theta = (0.1)^2;
        acr1 = 0.001;             
        display(' ');        
        display(['substantial heterogeneity case: theta=[',num2str(thetalb),',',num2str(thetaub),'], psi=',num2str(psi),', phi=',num2str(phi),', zeta=',num2str(zeta),'.']);        
    end  

    theta = thetalb:acr1:thetaub; 
    [Theta,W] = meshgrid(theta,w); 
    X = [Theta(:) W(:)];      
    mu = [mu_theta mu_w];      
    Sigma = [eps2_theta psi*(eps2_theta^0.5)*(eps2_w^0.5); psi*(eps2_theta^0.5)*(eps2_w^0.5) eps2_w];
    y = mvnpdf(X,mu,Sigma); % the CDF in Tx1
    y = reshape(y,length(w),length(theta)); % change to NxM matrix
    pdf = y; 
    tp = 1- sum(sum(y))*acr1*acr2; % trancated percentage
    Etheta = sum(sum(Theta.*y))*acr1*acr2; % check mean of theta
    Ew = sum(sum(W.*y))*acr1*acr2;    % check the mean of w
    thetawj = [mu_theta thetaub];  % to check if there are multiple equilibra; by continuety, the lower- and upper-end points are enough 
    
    g=0; % smaller g   
    for j=1:length(thetawj)
        thetaw = thetawj(j); % set a starting value for thetaw
        % display(thetaw); 
        idx1=0;  % idx1 controls for the first interation
        while (idx1==0) || (abs(thetawjj-thetaw) >= 0.00001)
            if idx1 == 1
                thetaw = thetawjj;
            end 
            Ai = (1+r)/(1+r+thetaw+g-thetaw*g);
            c1 = W./(1+thetaw*(((1+r)*Theta/(1+rho)/thetaw).^(1/phi))/(1+r) + ...
                thetaw*(((1+r)*Theta/(1+rho)/thetaw).^(1/phi))*(((1+r)*zeta/(1+rho))^(1/phi))/(1+r)/(1+r) ...
                + (1-thetaw)*(((1+r)*(1-Theta).*zeta/(1+rho)/(1-thetaw)).^(1/phi))/(1+r)); 
            c2 = (((1+r)*Theta/(1+rho)/thetaw).^(1/phi)).*c1;
            b2 = (((1+r)*(1-Theta).*zeta/(1+rho)/(1-thetaw)).^(1/phi)).*c1;     
            b3 = ((((1+r)*Theta/(1+rho)/thetaw).^(1/phi))*(((1+r)*zeta/(1+rho))^(1/phi))).*c1;
            alpha = (c2+b3/(1+r)-b2)./((1-g)*Ai);     
            alpha(alpha(:)<0)=0; 
            alpha(alpha(:)>m)=m;

            fun_n= Theta.*alpha.*y;
            fun_d= alpha.*y;
            thetawjj = sum(sum(fun_n)) / sum(sum(fun_d));     
            idx1 = 1; 
            % display(thetawjj)
        end
        if j==1
            thetaw00 = thetaw;  % to store the first equilibrium value if multiple equilibria exist or the only equilibrium value;     
        end 

        if (j>1) && abs(thetaw-thetaw00)>0.0001;   % if another equilibrium value thetaw exists 
            thetaw00m = thetaw;  % to store another equilibrium value
            display(['When g=',num2str(g),', multiple equilibria with' ,' thetaw=',num2str(thetaw00m), '! Ignore it and select the one with lowest thetaw instead.']);  
            if thetaw < thetaw00 
                thetaw00 = thetaw;  % Select the one with lower thetaw
            end
            pause(2)
            % stop % break
        end     
    end      
    alphag00 = alpha; A00 = Ai; 
    pause(2)
    display(['When g=',num2str(g),', thetaw=',num2str(thetaw00),'.']);      


    g=0.20;     % larger g
    for j=1:length(thetawj)
        thetaw = thetawj(j); % set a starting value for thetaw
        % display(thetaw); 
        idx1=0;  % idx1 controls for the first interation
        while (idx1==0) || (abs(thetawjj-thetaw) >= 0.00001)
            if idx1 == 1
                thetaw = thetawjj;
            end 
            Ai = (1+r)/(1+r+thetaw+g-thetaw*g);
            c1 = W./(1+thetaw*(((1+r)*Theta/(1+rho)/thetaw).^(1/phi))/(1+r) + ...
                thetaw*(((1+r)*Theta/(1+rho)/thetaw).^(1/phi))*(((1+r)*zeta/(1+rho))^(1/phi))/(1+r)/(1+r) ...
                + (1-thetaw)*(((1+r)*(1-Theta).*zeta/(1+rho)/(1-thetaw)).^(1/phi))/(1+r)); 
            c2 = (((1+r)*Theta/(1+rho)/thetaw).^(1/phi)).*c1;
            b2 = (((1+r)*(1-Theta).*zeta/(1+rho)/(1-thetaw)).^(1/phi)).*c1;     
            b3 = ((((1+r)*Theta/(1+rho)/thetaw).^(1/phi))*(((1+r)*zeta/(1+rho))^(1/phi))).*c1;
            alpha = (c2+b3/(1+r)-b2)./((1-g)*Ai);     
            alpha(alpha(:)<0)=0; 
            alpha(alpha(:)>m)=m;

            fun_n= Theta.*alpha.*y;
            fun_d= alpha.*y;
            thetawjj = sum(sum(fun_n)) / sum(sum(fun_d));     
            idx1 = 1; 
        end
        if j==1
            thetaw20 = thetaw;        
        end 

        if (j>1) && abs(thetaw-thetaw20)>0.0001;  
            thetaw20m = thetaw;  % to store another equilibrium value            
            display(['When g=',num2str(g),', multiple equilibria with',' thetaw=',num2str(thetaw20m), '! Ignore it and select the one with lowest thetaw instead.']);  
            if thetaw < thetaw20 
                thetaw20 = thetaw;  % Select the one with lower thetaw
            end            
            pause(2)
            % stop % break
        end     
    end
    alphag20 = alpha; A20 = Ai; 
    display(['When g=',num2str(g),', thetaw=',num2str(thetaw20),'.']);          
   
%    display([tp A00 A20 thetaw00 thetaw20 thetaw00-thetaw20]);
    diff = round((thetaw20-thetaw00)*10000)/10000;
    display(['The difference is ', num2str(diff),'.']);   
    pause(2)    
    
    if i == 1 
        thetaA = theta; alphag00A= alphag00; alphag20A = alphag20; A00A = A00; A20A = A20; 
    elseif i == 2
        thetaC = theta; alphag00C= alphag00; alphag20C = alphag20; A00C = A00; A20C = A20;             
    end    
    
end    



